Hunting for TFs using tradeSeq: incorporating derivative threshold
Koen Van den Berge
1/17/2020
here::set_here()## File .here already exists in /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/scriptssuppressPackageStartupMessages({
library(slingshot)
library(tradeSeq)
library(SingleCellExperiment)
library(cowplot)
library(rgl)
library(clusterExperiment)
library(RColorBrewer)
library(aggregation)
library(ggplot2)
library(pheatmap)
library(wesanderson)
library(UpSetR)
library(gridExtra)
})Import data
sds <- readRDS("../data/finalTrajectory/sling.rds")
counts <- readRDS("../data/finalTrajectory/counts_noResp_noMV.rds")
counts <- round(counts)
sce <- readRDS("../data/finalTrajectory/sce_tradeSeq20200904.rds")
load("../data/ALL_TF.Rda")
clDatta <- readRDS("../data/finalTrajectory/dattaCl_noResp_noMV.rds")Helper functions
source("https://raw.githubusercontent.com/statOmics/tradeSeq/master/R/utils.R")
epsilon <- 1e-6
calculateDerivatives <- function(sce, gene, lineage, epsilon=1e-6, nPoints=100){
## gather all required objects
id <- which(rownames(sce) == gene)
beta <- matrix(unlist(rowData(sce)$tradeSeq$beta[id,]),ncol=1)
if(any(is.na(beta))){
dfsd <- data.frame(der=rep(NA, nPoints), sdDer=rep(NA, nPoints))
colnames(dfsd) <- c(paste0("der",lineage), paste0("sdDer",lineage))
return(dfsd)
}
Sigma <- rowData(sce)$tradeSeq$Sigma[[id]]
X <- colData(sce)$tradeSeq$X
slingshotColData <- colData(sce)$slingshot
pseudotime <- slingshotColData[,grep(x = colnames(slingshotColData),
pattern = "pseudotime")]
dm <- colData(sce)$tradeSeq$dm
## calculate lp matrix with finite differencing
newd1 <- .getPredictRangeDf(dm, lineage, nPoints=nPoints)
newd2 <- newd1
newd1[,paste0("t",lineage)] <- newd1[,paste0("t",lineage)] - epsilon
X0 <- predictGAM(X, newd1, pseudotime)
newd2[,paste0("t",lineage)] <- newd2[,paste0("t",lineage)] + epsilon
X1 <- predictGAM(X, newd2, pseudotime)
Xp <- (X1-X0)/(2*epsilon) ## maps coefficients to (fd approx.) derivatives
## calculate first derivative
df <- Xp%*%beta ## ith smooth derivative
df.sd <- rowSums(Xp%*%Sigma*Xp)^.5 ## cheap diag(Xi%*%b$Vp%*%t(Xi))^.5
dfsd <- data.frame(der=df, sdDer=df.sd)
colnames(dfsd) <- c(paste0("der",lineage), paste0("sdDer",lineage))
return(dfsd)
}
calculateDerivativesWithSigma <- function(sce, gene, lineage, epsilon=1e-6, nPoints=100){
## gather all required objects
id <- which(rownames(sce) == gene)
beta <- matrix(unlist(rowData(sce)$tradeSeq$beta[id,]),ncol=1)
Sigma <- rowData(sce)$tradeSeq$Sigma[[id]]
X <- colData(sce)$tradeSeq$X
slingshotColData <- colData(sce)$slingshot
pseudotime <- slingshotColData[,grep(x = colnames(slingshotColData),
pattern = "pseudotime")]
dm <- colData(sce)$tradeSeq$dm
## calculate lp matrix with finite differencing
newd1 <- .getPredictRangeDf(dm, lineage, nPoints=nPoints)
newd2 <- newd1
newd1[,paste0("t",lineage)] <- newd1[,paste0("t",lineage)] - epsilon
X0 <- predictGAM(X, newd1, pseudotime)
newd2[,paste0("t",lineage)] <- newd2[,paste0("t",lineage)] + epsilon
X1 <- predictGAM(X, newd2, pseudotime)
Xp <- (X1-X0)/(2*epsilon) ## maps coefficients to (fd approx.) derivatives
## calculate first derivative
df <- Xp%*%beta ## ith smooth derivative
#df.sd <- rowSums(Xp%*%Sigma*Xp)^.5 ## cheap diag(Xi%*%b$Vp%*%t(Xi))^.5
sigmaDer <- Xp%*%Sigma%*%t(Xp)
dfsd <- list(der=df, sdDer=sigmaDer)
names(dfsd) <- c(paste0("der",lineage), paste0("sdDer",lineage))
return(dfsd)
}
plotSmoothersWithDerivatives <- function(sce, gene, counts, gridDf, derDf){
dfAll <- as.data.frame(cbind(gridDf,derDf))
cowplot::plot_grid(
#smoothers
plotSmoothers(sce, gene=gene, counts=counts, lwd=1),
ggplot(dfAll, aes(x=grid1, y=derDf$der1)) +
geom_line() +
geom_line(aes(x=grid1, y=derDf$der1 - 2*derDf$sdDer1), lty=2) +
geom_line(aes(x=grid1, y=derDf$der1 + 2*derDf$sdDer1), lty=2) +
geom_hline(yintercept=0, col="red") +
ggtitle("Neuronal lineage") +
ylab("First derivative"),
ggplot(dfAll, aes(x=grid2, y=derDf$der2)) +
geom_line() +
geom_line(aes(x=grid2, y=derDf$der2 - 2*derDf$sdDer2), lty=2) +
geom_line(aes(x=grid2, y=derDf$der2 + 2*derDf$sdDer2), lty=2) +
geom_hline(yintercept=0, col="red") +
ggtitle("rHBC lineage") +
ylab("First derivative"),
ggplot(dfAll, aes(x=grid3, y=derDf$der3)) +
geom_line() +
geom_line(aes(x=grid3, y=derDf$der3 - 2*derDf$sdDer3), lty=2) +
geom_line(aes(x=grid3, y=derDf$der3 + 2*derDf$sdDer3), lty=2) +
geom_hline(yintercept=0, col="red") +
ggtitle("SUS lineage") +
ylab("First derivative"),
# test statistics
ggplot(dfAll, aes(x=grid1, y=derDf$der1/derDf$sdDer1)) +
geom_line() +
geom_line(aes(x=grid2, y=derDf$der2/derDf$sdDer2)) +
geom_line(aes(x=grid3, y=derDf$der3/derDf$sdDer3)) +
geom_hline(yintercept=0, col="red") +
ggtitle("Test statistics"),
nrow=3, ncol=2
)
}
makeGSEATable <- function(overlapFile, nBP=20){
overlap <- readLines(overlapFile)
overlapSets <- overlap[10:30]
overlapSetsSplit <- sapply(overlapSets, function(x) strsplit(x, split = "\t"))
gsNames <- unname(unlist(lapply(overlapSetsSplit, "[[", 1)))
genesInSet <- unname(unlist(lapply(overlapSetsSplit, "[[", 2)))
genesInOverlap <- unname(unlist(lapply(overlapSetsSplit, "[[", 4)))
pvalUniqGam <- unname(unlist(lapply(overlapSetsSplit, "[[", 6)))
qval <- unname(unlist(lapply(overlapSetsSplit, "[[", 7)))
gsNames <- tolower(gsNames)
gsNames <- gsub(x = gsNames, pattern = "_", replacement = " ")
gsNames[-1] <- unname(sapply(gsNames[-1], function(x) substr(x, 4, nchar(x))))
tab <- data.frame(
geneSet = gsNames[-1],
overlap = as.numeric(genesInOverlap[-1]),
genesInSet = as.numeric(genesInSet[-1]),
qvalue = qval[-1]
)
#library(xtable)
#print(xtable(tab), type="html")
return(tab)
}Estimate first derivatives
tf <- intersect(ALL_TF, rownames(sce))
# derAllTF <- list()
# for(ii in 1:length(tf)){
# gene <- tf[ii]
# derGene <- do.call(cbind, sapply(1:3, function(ll){
# calculateDerivatives(sce=sce,
# gene=gene,
# lineage=ll,
# nPoints=100)
# }, simplify=FALSE))
# derAllTF[[ii]] <- derGene
# }
# saveRDS(derAllTF, file="../data/derAllTF_20200916.rds")
derAllTF <- readRDS("../data/derAllTF_20200916.rds")
names(derAllTF) <- tfOrder TFs based on expression peak, require a threshold on derivative to be called significant
A gene’s expression peaks when its first derivative equals zero. Hence, we could use this to order the TFs which have a significant expression peak.
Since the goal of this analysis is ordering of genes, we will ignore multiple testing correction across genes. However, since we want to restrict ourselves to genes that actually do have a significant expression peak, we will do within-gene FWER correction.
nPoints <- 100
dm <- colData(sce)$tradeSeq$dm
grid1 <- tradeSeq:::.getPredictRangeDf(dm = dm,
lineageId = 1,
conditionId = NULL,
nPoints = nPoints)
cp <- c('#1B9E77', '#E6AB02', '#E7298A', '#66A61E', '#BEAED4' ,'#D95F02', '#A6761D', '#666666', '#1F78B4') #Rebecca colors
cl <- factor(colnames(sds@clusterLabels)[apply(sds@clusterLabels,1,which.max)])
colDf <- data.frame(cl=levels(cl),
rcCol=cp, #RKC colors
trCol=brewer.pal(9,'Set1')) #KS colors
### write a function that creates an Anno object with length as nPoints.
anno <- function(lineage, nPoints, sds, dm, cl, lengthThresh=5){
cp <- c('#1B9E77', '#E6AB02', '#E7298A', '#66A61E', '#BEAED4' ,'#D95F02', '#A6761D', '#666666', '#1F78B4') #Rebecca colors
#cl <- factor(colnames(sds@clusterLabels)[apply(sds@clusterLabels,1,which.max)])
pt <- slingPseudotime(sds)[,lineage]
#cw <- slingCurveWeights(sds)
#linID <- apply(cw,1,which.max)
linID <- apply(dm[,paste0("l",1:3)], 1, which.max)
pt1 <- pt[linID == lineage]
cl1 <- droplevels(cl[linID == lineage])
## divide pseudotime in nPoint
#qq <- quantile(pt1, probs = seq(0,1,length=nPoints+1))
qq <- seq(0, max(pt1), length.out=nPoints+1)
allCols <- c()
for(kk in 1:nPoints){
id <- which(pt1 > qq[kk] & pt1 < qq[kk+1])
if(length(id) < lengthThresh){
allCols[kk] <- "white"
next
}
maxCl <- names(sort(table(droplevels(cl1[id])), decreasing=TRUE))[1]
col <- which(levels(cl) == maxCl)
allCols[kk] <- cp[col]
}
return(allCols)
}
orderPeakGenesThresh <- function(grid1, derivatives, tf, lineage, sce, cl,
main=NULL, ablines=NULL, plotHist=TRUE,
plotLines = TRUE, plotHeatmap = TRUE,
threshold = 0.025, clusteredHeatmap=TRUE,
show_rownames = TRUE, show_colnames=FALSE,
returnOnlyHeatmap = FALSE){
# grid1 is the grid on which the derivatives have been calculated
# derivatives is a list of first derivatives and their SD
# tf is a vector of all TF names
# lineage is the lineage of interest
# sce are the fitted tradeSeq models
# cl are the cluster labels
nPoints <- nrow(grid1)
# 1. first select genes with a significant first derivative anywhere in the lineage.
# I do not consider negative first derivatives here since we're focussed on peaks.
## a. calculate p-values without threshold => for selecting uniquely peaking genes
pvalNeurAllNoThresh <- lapply(derivatives, function(x){
der <- x[,(lineage*2 - 1)]
derSD <- x[,(lineage*2)]
testStat <- der / derSD
# testStat <- (abs(der)-threshold) / derSD
pvals <- pmin(1,(1-pnorm(testStat)))
return(pvals)
})
pvalNoThresh <- do.call(rbind,pvalNeurAllNoThresh)
## b. calculate p-values with threshold.
resNeurAll <- lapply(derivatives, function(x){
der <- x[,(lineage*2 - 1)]
derSD <- x[,(lineage*2)]
# testStat <- der / derSD
testStat <- (abs(der)-threshold) / derSD
pvals <- pmin(1,(1-pnorm(testStat)))
pvals[der < threshold] <- 1
return(list(testStat=testStat,
pvals=pvals))
})
pvalNeurAll <- lapply(resNeurAll, function(x) x$pvals)
testNeurAll <- lapply(resNeurAll, function(x) x$testStat)
pvalThresh <- do.call(rbind,pvalNeurAll)
testThresh <- do.call(rbind,testNeurAll)
maxTestsThresh <- matrixStats::rowMaxs(testThresh)
# Check significance.
neurPeakGenes <- tf[unlist(lapply(pvalNeurAll, function(x) any(p.adjust(x, "holm") <= 0.05)))]
pvalNeurAll <- pvalNeurAll[neurPeakGenes]
length(neurPeakGenes) / length(tf)
# 2. then define the most pronounced peak by looking at the maximum test statistic for each gene
tLin <- grid1[,paste0("t",lineage)]
tNeur <- tLin[unlist(lapply(pvalNeurAll, which.min))]
# hist(tNeur)
# 3. then define where that peak peaks by checking where the first derivative crosses zero the first time AFTER the pseudotime value defined above.
firstPeak <- c()
for(ii in 1:length(neurPeakGenes)){
gene <- neurPeakGenes[ii]
tMax <- tNeur[ii]
d1 <- derivatives[[gene]][,paste0("der",lineage)]
# check where the derivative crosses zero first time AFTER tNeur
d1 <- d1[grid1[,paste0("t",lineage)] > tMax]
t1 <- grid1[,paste0("t",lineage)][grid1[,paste0("t",lineage)] > tMax]
tCrossZero <- t1[which(diff(sign(d1))==-2)]
if(length(tCrossZero)==0){
firstPeak[ii] <- NA
next
}
firstPeak[ii] <- min(tCrossZero)
}
names(firstPeak) <- neurPeakGenes
firstPeak[is.na(firstPeak)] <- max(grid1)
oo <- order(firstPeak, decreasing=FALSE)
if(plotHist){
hist(firstPeak, breaks=50, main=paste0("Peak times: ", main),
xlim = c(0,max(grid1)))
}
## line plot
X <- colData(sce)$tradeSeq$X # linear predictor
slingshotColData <- colData(sce)$slingshot
pseudotime <- slingshotColData[,grep(x = colnames(slingshotColData),
pattern = "pseudotime")]
betaMat <- rowData(sce)$tradeSeq$beta[[1]]
rownames(betaMat) <- rownames(sce)
Xdf <- predictGAM(lpmatrix = X,
df = grid1,
pseudotime = pseudotime)
yhatMat <- matrix(NA, nrow=length(firstPeak), ncol=nPoints)
rownames(yhatMat) <- names(firstPeak)
for(ii in 1:length(firstPeak)){
beta <- betaMat[rownames(yhatMat)[ii],]
yhat <- c(exp(t(Xdf %*% t(beta))))
yhatMat[ii,] <- yhat
}
yhatMatScaled <- t(scale(t(yhatMat)))
if(plotLines){
pal <- wesanderson::wes_palette("Zissou1", n=length(firstPeak), type="continuous")
plot(x=grid1[,paste0("t",lineage)], y=seq(min(yhatMatScaled), max(yhatMatScaled), length.out=nPoints), type='n', bty='l',
ylab="Standardized Expression", xlab="Pseudotime", main=main)
for(ii in 1:length(firstPeak)){
lines(x=grid1[,paste0("t",lineage)], y=yhatMatScaled[names(firstPeak)[oo][ii],], col=pal[ii])
}
if(!is.null(ablines)) abline(v=ablines, col="black", lwd=3)
}
## heatmap
if(plotHeatmap){
annCols <- anno(lineage, nPoints, sds, dm, cl)
annCols2 <- as.character(colDf$trCol[match(annCols, colDf$rcCol)])
annCols2[is.na(annCols2)] <- "white"
dfAnn <- data.frame(cellType=annCols, cellType2=annCols2)
rownames(dfAnn) <- colnames(yhatMatScaled)
annColors <- list()
annColors[[1]] <- as.character(levels(factor(annCols)))
names(annColors[[1]]) <- levels(factor(annCols))
annColors[[2]] <- as.character(levels(factor(annCols2)))
names(annColors[[2]]) <- levels(factor(annCols2))
names(annColors) <- c("cellType", "cellType2")
pal <- wesanderson::wes_palette("Zissou1", n=12, type="continuous")
yhatMatScaledOrdered <- yhatMatScaled[names(firstPeak)[oo],]
colnames(yhatMatScaledOrdered) <- rownames(dfAnn) <- as.character(1:ncol(yhatMatScaledOrdered))
ph <- pheatmap(yhatMatScaledOrdered, cluster_cols=FALSE, cluster_rows=FALSE,
border_color=NA, col=pal, main=main, annotation_col=dfAnn[,2,drop=FALSE],
annotation_colors = annColors, annotation_names_col = FALSE,
annotation_legend = FALSE, show_rownames = show_rownames,
show_colnames = show_colnames)#, colors=cols, breaks=breaks)
ph
if(clusteredHeatmap){
ph <- pheatmap(yhatMatScaledOrdered, cluster_cols=FALSE, cluster_rows=TRUE,
border_color=NA, col=pal, main=paste0(main,": Clustered genes"),
annotation_col=dfAnn[,2,drop=FALSE],
annotation_colors = annColors, annotation_names_col = FALSE,
annotation_legend = FALSE, show_rownames = show_rownames,
show_colnames = show_colnames)
ph
}
}
if(returnOnlyHeatmap){
return(ph)
} else {
return(list(firstPeak=firstPeak, standExpr=yhatMatScaled,
yhat=yhatMat,
pvalNoThresh=pvalNoThresh, pvalThresh=pvalThresh,
maxTestStats = maxTestsThresh))
}
}
heatmapScaledAcrossAllLineages <- function(models, genes, sds, nPoints, outFile, g=20,
cl, showRowNames=TRUE, width=40, height=20,
showLegend = TRUE){
### scaled across all lineages
yhat <- tradeSeq::predictSmooth(models=models, gene=genes, nPoints=nPoints, tidy=FALSE)
yhat_stand <- t(scale(t(log(yhat))))
titles <- c("Neuronal", "Sustentacular", "rHBC")
yhatList <- list(neur=yhat_stand[,1:nPoints],
rhbc=yhat_stand[,(nPoints+1):(2*nPoints)],
sus=yhat_stand[,(2*nPoints+1):(3*nPoints)])
dm <- colData(models)$tradeSeq$dm
pal <- wesanderson::wes_palette("Zissou1", n=g+1, type="continuous")
breaks <- Hmisc::cut2(yhat_stand, g=g+1, onlycuts=TRUE)
plotList <- list()
for(ii in 1:3){
annCols <- anno(ii, nPoints, sds, dm, cl)
annCols2 <- as.character(colDf$trCol[match(annCols, colDf$rcCol)])
annCols2[is.na(annCols2)] <- "white"
dfAnn <- data.frame(cellType=annCols, cellType2=annCols2)
rownames(dfAnn) <- colnames(yhatList[[ii]])
annColors <- list()
annColors[[1]] <- as.character(levels(factor(annCols)))
names(annColors[[1]]) <- levels(factor(annCols))
annColors[[2]] <- as.character(levels(factor(annCols2)))
names(annColors[[2]]) <- levels(factor(annCols2))
names(annColors) <- c("cellType", "cellType2")
ph <- pheatmap(yhatList[[ii]], cluster_cols=FALSE, cluster_rows=FALSE,
border_color=NA, main=titles[ii], breaks=breaks, color=pal,
show_colnames = FALSE, annotation_col=dfAnn[,2,drop=FALSE],
annotation_colors = annColors, annotation_names_col = FALSE,
annotation_legend = FALSE, show_rownames = showRowNames,
legend = showLegend)
plotList[[ii]] <- ph[[4]]
}
pp <- cowplot::plot_grid(plotlist=plotList, nrow=1)
ggsave(filename = outFile, plot = pp,
units = "in", width=width, height=height)
}Example plot of derivatives and associated test statistics
grid1 <- tradeSeq:::.getPredictRangeDf(dm = dm,
lineageId = 1,
conditionId = NULL,
nPoints = nPoints)
grid2 <- tradeSeq:::.getPredictRangeDf(dm = dm,
lineageId = 2,
conditionId = NULL,
nPoints = nPoints)
grid3 <- tradeSeq:::.getPredictRangeDf(dm = dm,
lineageId = 3,
conditionId = NULL,
nPoints = nPoints)
gridDf <- data.frame(grid1=grid1$t1,
grid2=grid2$t2,
grid3=grid3$t3)
plotSmoothersWithDerivatives(sce, "Junb", counts, gridDf = gridDf, derDf=derAllTF[["Junb"]])
Lineage-specific ordering
# neuronal lineage
pdf("../plots/peakAnalysis_Thresholded/neuronal/orderedPeaksPlots.pdf")
neurRes <- orderPeakGenesThresh(grid=grid1, derivatives=derAllTF, tf=tf, lineage=1, sce=sce,
cl=clDatta, main="Neuronal", threshold = 0.1, show_rownames = FALSE)
neurPeaks <- neurRes$firstPeak
dev.off()## pdf
## 3pdf("../plots/peakAnalysis_Thresholded/neuronal/orderedPeaksPlots_tallHeatmaps.pdf", width=10, height=65)
neurRes <- orderPeakGenesThresh(grid=grid1, derivatives=derAllTF, tf=tf, lineage=1, sce=sce,
cl=clDatta, main="Neuronal", plotHist=FALSE, plotLines=FALSE,
threshold = 0.1)
dev.off()## pdf
## 3# Sus lineage
pdf("../plots/peakAnalysis_Thresholded/sus/orderedPeaksPlots.pdf")
susRes <- orderPeakGenesThresh(grid1=grid2, derivatives=derAllTF, tf=tf, lineage=2, sce=sce,
cl=clDatta, main="Sustentacular", threshold = 0.1)
susPeaks <- susRes$firstPeak
dev.off()## pdf
## 3pdf("../plots/peakAnalysis_Thresholded/sus/orderedPeaksPlots_tallHeatmaps.pdf", width=10, height=40)
susRes <- orderPeakGenesThresh(grid1=grid2, derivatives=derAllTF, tf=tf, lineage=2, sce=sce,
cl=clDatta, main="Sustentacular", plotHist=FALSE, plotLines=FALSE, threshold = 0.1)
dev.off()## pdf
## 3# rHBC lineage
pdf("../plots/peakAnalysis_Thresholded/sus/orderedPeaksPlots.pdf")
rhbcRes <- orderPeakGenesThresh(grid1=grid3, derivatives=derAllTF, tf=tf, lineage=3, sce=sce,
cl=clDatta, main="rHBC", threshold = 0.1)
rhbcPeaks <- rhbcRes$firstPeak
dev.off()## pdf
## 3pdf("../plots/peakAnalysis_Thresholded/rhbc/orderedPeaksPlots_tallHeatmaps.pdf", width=10, height=60)
rhbcRes <- orderPeakGenesThresh(grid1=grid3, derivatives=derAllTF, tf=tf, lineage=3, sce=sce,
cl=clDatta, main="rHBC", plotHist=FALSE, plotLines=FALSE,
threshold = 0.1)
dev.off()## pdf
## 3peakList <- list(neur=neurPeaks,
sus=susPeaks,
rhbc=rhbcPeaks)
peakResList <- list(neur=neurRes,
sus=susRes,
rhbc=rhbcRes)
lapply(peakList, length)## $neur
## [1] 524
##
## $sus
## [1] 231
##
## $rhbc
## [1] 284saveRDS(peakList, file="../data/peakList_thresholded.rds")
saveRDS(peakResList, file="../data/peakRes_thresholded.rds")Custom plots for paper
# Neuronal
png("../plots/peakAnalysis_Thresholded/neuronal/lines_neuronal.png", width=9, height=10, units="in", res=300)
neurRes <- orderPeakGenesThresh(grid=grid1, derivatives=derAllTF, tf=tf, lineage=1, sce=sce,
cl=clDatta, main="Neuronal", plotHist=FALSE, plotLines=TRUE,
plotHeatmap = FALSE, threshold = 0.1, clusteredHeatmap = FALSE,
show_rownames=FALSE)
dev.off()## quartz_off_screen
## 2png("../plots/peakAnalysis_Thresholded/neuronal/cascadeHeatmap_neuronal.png", width=9, height=10, units="in", res=300)
neurRes <- orderPeakGenesThresh(grid=grid1, derivatives=derAllTF, tf=tf, lineage=1, sce=sce,
cl=clDatta, main="Neuronal", plotHist=FALSE, plotLines=FALSE,
plotHeatmap = TRUE, threshold = 0.1, clusteredHeatmap = FALSE,
show_rownames=FALSE)
dev.off()## pdf
## 3# Sus
png("../plots/peakAnalysis_Thresholded/sus/lines_sus.png", width=9, height=10, units="in", res=300)
susRes <- orderPeakGenesThresh(grid=grid2, derivatives=derAllTF, tf=tf, lineage=2, sce=sce,
cl=clDatta, main="Sustentacular", plotHist=FALSE, plotLines=TRUE,
plotHeatmap = FALSE, threshold = 0.1, clusteredHeatmap = FALSE,
show_rownames=FALSE)
dev.off()## pdf
## 3png("../plots/peakAnalysis_Thresholded/sus/cascadeHeatmap_sus.png", width=9, height=10, units="in", res=300)
susRes <- orderPeakGenesThresh(grid=grid2, derivatives=derAllTF, tf=tf, lineage=2, sce=sce,
cl=clDatta, main="Sustentacular", plotHist=FALSE, plotLines=FALSE,
plotHeatmap = TRUE, threshold = 0.1, clusteredHeatmap = FALSE,
show_rownames=FALSE)
dev.off()## pdf
## 3# HBC
png("../plots/peakAnalysis_Thresholded/rhbc/lines_rhbc.png", width=9, height=10, units="in", res=300)
rhbcRes <- orderPeakGenesThresh(grid=grid3, derivatives=derAllTF, tf=tf, lineage=3, sce=sce,
cl=clDatta, main="rHBC", plotHist=FALSE, plotLines=TRUE,
plotHeatmap = FALSE, threshold = 0.1, clusteredHeatmap = FALSE,
show_rownames=FALSE)
dev.off()## pdf
## 3png("../plots/peakAnalysis_Thresholded/rhbc/cascadeHeatmap_rhbc.png", width=9, height=10, units="in", res=300)
rhbcRes <- orderPeakGenesThresh(grid=grid3, derivatives=derAllTF, tf=tf, lineage=3, sce=sce,
cl=clDatta, main="rHBC", plotHist=FALSE, plotLines=FALSE,
plotHeatmap = TRUE, threshold = 0.1, clusteredHeatmap = FALSE,
show_rownames=FALSE)
dev.off()## pdf
## 3All cascades next to each other
phNeur <- orderPeakGenesThresh(grid=grid1, derivatives=derAllTF, tf=tf, lineage=1, sce=sce,
cl=clDatta, main="Neuronal", plotHist=FALSE, plotLines=FALSE,
plotHeatmap = TRUE, threshold = 0.1, clusteredHeatmap = FALSE,
show_rownames=FALSE, returnOnlyHeatmap = TRUE)
phSus <- orderPeakGenesThresh(grid=grid2, derivatives=derAllTF, tf=tf, lineage=2, sce=sce,
cl=clDatta, main="Sustentacular", plotHist=FALSE, plotLines=FALSE,
plotHeatmap = TRUE, threshold = 0.1, clusteredHeatmap = FALSE,
show_rownames=FALSE, returnOnlyHeatmap = TRUE)
phHBC <- orderPeakGenesThresh(grid=grid3, derivatives=derAllTF, tf=tf, lineage=3, sce=sce,
cl=clDatta, main="rHBC", plotHist=FALSE, plotLines=FALSE,
plotHeatmap = TRUE, threshold = 0.1, clusteredHeatmap = FALSE,
show_rownames=FALSE, returnOnlyHeatmap = TRUE)
cowplot::plot_grid(phNeur[[4]], phSus[[4]], phHBC[[4]],
nrow=1, ncol=3)
ggsave("../plots/fig2_allCascades.png")## Saving 7 x 5 in imageIdentify shared and unique lineage TFs
We can use the results above to identify shared and unique TFs for each lineage, using a set diff. For shared TFs, we first aggregate the p-values testing whether a derivative is different from zero within the gene, using Sidak aggregation. Next, we perform FDR correction on the aggregated p-values for each lineage separately. The TFs that are significant post FDR correction in all three lineages are considered to be shared.
For TFs uniquely peaking in a lineage, we retain TFs that are significantly peaking in that lineage, and that are higher expressed in that lineage versus the other two lineages, by a fold change threshold of \(1.5\) on the maximum fitted expression value.
yhatAllTF <- predictSmooth(sce, gene=tf, nPoints=100, tidy=FALSE)
library(aggregation)
## shared TFs
# aggregate p-values using Sidak method for each lineage
aggPvals <- cbind(apply(neurRes$pvalThresh, 1, sidak),
apply(rhbcRes$pvalThresh, 1, sidak),
apply(susRes$pvalThresh, 1, sidak))## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value
## Warning in FUN(newX[, i], ...): Extreme p-values around and below 10e-17 will
## produce aggregated p-value of 0. Replace extreme p-values with 10e-17 to obtain
## upper bound on the aggregated p-value# FDR correction for each lineage across all TFs
aggPvalsFDR <- apply(aggPvals,2,p.adjust,method="fdr")
barplot(table(rowSums(aggPvalsFDR <= 0.05)), xlab="Number of lineages in which a TF significantly peaks.")
sharedTF <- rownames(aggPvalsFDR)[which(rowSums(aggPvalsFDR <= 0.05) == 3)]
length(sharedTF)## [1] 65## unique TFs
getUniqueTFs <- function(peakResList, yhatAllTF, lineage, name, fcThreshold=1.5){
# get TFs identified for that lineage
tfLineage <- names(peakResList[[name]]$firstPeak)
# within that set, check if expression is much higher
linID <- ((lineage-1)*100+1):(lineage*100)
delta <- rowMax(yhatAllTF[tfLineage,linID]) / rowMax(yhatAllTF[tfLineage,-linID])
tfLineageDelta <- tfLineage[delta > fcThreshold]
return(tfLineageDelta)
}
neurUniqTFs <- getUniqueTFs(peakResList, yhatAllTF, lineage=1, name="neur")
susUniqTFs <- getUniqueTFs(peakResList, yhatAllTF, lineage=2, name="sus")
rhbcUniqTFs <- getUniqueTFs(peakResList, yhatAllTF, lineage=3, name="rhbc")
length(neurUniqTFs) ; length(susUniqTFs) ; length(rhbcUniqTFs) ## [1] 352## [1] 31## [1] 61# write.table(sharedTF, file="../data/sharedTF.txt", row.names=FALSE, col.names=FALSE, quote=FALSE)
# write.table(neurUniqTFs, file="../data/neurUniqTf_threshold.txt", row.names=FALSE, col.names=FALSE, quote=FALSE)
# write.table(susUniqTFs, file="../data/susUniqTf_threshold.txt", row.names=FALSE, col.names=FALSE, quote=FALSE)
# write.table(rhbcUniqTFs, file="../data/rhbcUniqTf_threshold.txt", row.names=FALSE, col.names=FALSE, quote=FALSE)Alternative approach to identify shared TFs using Fisher p-value aggregation
yhatAllTF <- predictSmooth(sce, gene=tf, nPoints=100, tidy=FALSE)
fisherX <- function(pval){
pval[pval==0] <- 1e-100
stat <- -2 * sum(log(pval))
return(stat)
}
## shared TFs
# aggregate p-values using Fisher method for each lineage
aggStats <- cbind(apply(neurRes$pvalThresh, 1, fisherX),
apply(rhbcRes$pvalThresh, 1, fisherX),
apply(susRes$pvalThresh, 1, fisherX))
plot(density(log(aggStats[,1])), ylim=c(0,0.4))
lines(density(log(aggStats[,2])), col="blue")
lines(density(log(aggStats[,3])), col="darkseagreen3")
abline(v=7.5, col="red")
sharedTF1 <- rownames(aggStats)[which(rowSums(log(aggStats) > 7.5) == 3)]
sharedTF1## [1] "Ebf1" "Lhx2" "Fos" "Nupr1"## unique TFs
getUniqueTFs <- function(peakResList, yhatAllTF, lineage, name, fcThreshold=1.5){
# get TFs identified for that lineage
tfLineage <- names(peakResList[[name]]$firstPeak)
# within that set, check if expression is much higher
linID <- ((lineage-1)*100+1):(lineage*100)
delta <- rowMax(yhatAllTF[tfLineage,linID]) / rowMax(yhatAllTF[tfLineage,-linID])
tfLineageDelta <- tfLineage[delta > fcThreshold]
return(tfLineageDelta)
}
neurUniqTFs <- getUniqueTFs(peakResList, yhatAllTF, lineage=1, name="neur")
susUniqTFs <- getUniqueTFs(peakResList, yhatAllTF, lineage=2, name="sus")
rhbcUniqTFs <- getUniqueTFs(peakResList, yhatAllTF, lineage=3, name="rhbc")
length(neurUniqTFs) ; length(susUniqTFs) ; length(rhbcUniqTFs) ## [1] 352## [1] 31## [1] 61write.table(sharedTF1, file="../data/sharedTF_v2.txt", row.names=FALSE, col.names=FALSE, quote=FALSE)
write.table(neurUniqTFs, file="../data/neurUniqTf_threshold.txt", row.names=FALSE, col.names=FALSE, quote=FALSE)
write.table(susUniqTFs, file="../data/susUniqTf_threshold.txt", row.names=FALSE, col.names=FALSE, quote=FALSE)
write.table(rhbcUniqTFs, file="../data/rhbcUniqTf_threshold.txt", row.names=FALSE, col.names=FALSE, quote=FALSE)Another alternative approach, using pseudotime agreement
# 90 TFs in all three
sharedTF1 <- intersect(intersect(names(neurRes$firstPeak), names(susRes$firstPeak)), names(rhbcRes$firstPeak))
# check for similar pseudotime
ptSharedTF1 <- data.frame(neur=neurRes$firstPeak[sharedTF1],
sus=susRes$firstPeak[sharedTF1],
rhbc=rhbcRes$firstPeak[sharedTF1])
deltaTF1 <- rowDiffs(rowRanges(as.matrix(ptSharedTF1)))
hist(deltaTF1, breaks=50)
plot(density(log(deltaTF1)))
sharedTFDelta <- sharedTF1[deltaTF1[,1] < 1]
heatmapScaledAcrossAllLineages(models=sce,
genes=sharedTFDelta,
nPoints=100,
sds=sds,
cl=clDatta,
height=5,
width=5,
showRowNames = FALSE,
showLegend = FALSE,
outFile="../plots/peakAnalysis_Thresholded/shared/sharedTF_scaledAcrossAllLineages_delta.pdf")
heatmapScaledAcrossAllLineages(models=sce,
genes=sharedTFDelta,
nPoints=100,
sds=sds,
cl=clDatta,
height=7,
width=5,
showRowNames = TRUE,
showLegend = FALSE,
outFile="../plots/peakAnalysis_Thresholded/shared/sharedTF_scaledAcrossAllLineages_delta_tfNames.pdf")
Heatmaps of shared and unique transcription factors
Unique TFs for each lineage
Tall heatmaps for interpretation as well as smaller heatmaps for paper.
Smoother plots of unique genes
pdf("../plots/peakAnalysis_Thresholded/neuronal/smootherPlots_neurTFs_v2.pdf")
for(ii in 1:length(neurUniqTFs)){
print(plotSmoothers(sce, counts, sort(neurUniqTFs)[ii], pointCol = clDatta) +
ggtitle(sort(neurUniqTFs)[ii]))
}
dev.off()## quartz_off_screen
## 2pdf("../plots/peakAnalysis_Thresholded/sus/smootherPlots_susTFs_v2.pdf")
for(ii in 1:length(susUniqTFs)){
print(plotSmoothers(sce, counts, sort(susUniqTFs)[ii], pointCol = clDatta) +
ggtitle(sort(susUniqTFs)[ii]))
}
dev.off()## quartz_off_screen
## 2pdf("../plots/peakAnalysis_Thresholded/rhbc/smootherPlots_rhbcTFs_v2.pdf")
for(ii in 1:length(rhbcUniqTFs)){
print(plotSmoothers(sce, counts, sort(rhbcUniqTFs)[ii], pointCol = clDatta) +
ggtitle(sort(rhbcUniqTFs)[ii]))
}
dev.off()## quartz_off_screen
## 2pdf("../plots/peakAnalysis_Thresholded/shared/smootherPlots_sharedTF_testStat.pdf")
for(ii in 1:length(sharedTF1)){
print(plotSmoothers(sce, counts, sort(sharedTF1)[ii], pointCol = clDatta) +
ggtitle(sort(sharedTF1)[ii]))
}
dev.off()## quartz_off_screen
## 2pdf("../plots/peakAnalysis_Thresholded/shared/smootherPlots_sharedTF_delta.pdf")
for(ii in 1:length(sharedTFDelta)){
print(plotSmoothers(sce, counts, sort(sharedTFDelta)[ii], pointCol = clDatta) +
ggtitle(sort(sharedTFDelta)[ii]))
}
dev.off()## quartz_off_screen
## 2Session Info
sessionInfo()## R version 3.6.2 (2019-12-12)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS Mojave 10.14.5
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] parallel stats4 stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] gridExtra_2.3 UpSetR_1.4.0
## [3] wesanderson_0.3.6 pheatmap_1.0.12
## [5] ggplot2_3.3.2 aggregation_1.0.1
## [7] RColorBrewer_1.1-2 clusterExperiment_2.6.1
## [9] bigmemory_4.5.36 rgl_0.100.30
## [11] cowplot_1.0.0 SingleCellExperiment_1.8.0
## [13] SummarizedExperiment_1.16.1 DelayedArray_0.12.2
## [15] BiocParallel_1.20.1 matrixStats_0.56.0
## [17] Biobase_2.46.0 GenomicRanges_1.38.0
## [19] GenomeInfoDb_1.22.0 IRanges_2.20.2
## [21] S4Vectors_0.24.3 BiocGenerics_0.32.0
## [23] tradeSeq_1.3.24 slingshot_1.4.0
## [25] princurve_2.1.4
##
## loaded via a namespace (and not attached):
## [1] utf8_1.1.4 tidyselect_1.1.0 RSQLite_2.2.0
## [4] AnnotationDbi_1.48.0 htmlwidgets_1.5.1 grid_3.6.2
## [7] combinat_0.0-8 docopt_0.6.1 Rtsne_0.15
## [10] RNeXML_2.4.2 munsell_0.5.0 codetools_0.2-16
## [13] miniUI_0.1.1.1 withr_2.5.0 colorspace_1.4-1
## [16] fastICA_1.2-2 knitr_1.28 uuid_0.1-2
## [19] rstudioapi_0.11 zinbwave_1.11.6 NMF_0.21.0
## [22] labeling_0.3 slam_0.1-47 GenomeInfoDbData_1.2.2
## [25] bit64_0.9-7 farver_2.0.3 rhdf5_2.30.1
## [28] rprojroot_1.3-2 vctrs_0.3.8 generics_0.0.2
## [31] xfun_0.26 R6_2.4.1 doParallel_1.0.15
## [34] VGAM_1.1-2 locfit_1.5-9.1 manipulateWidget_0.10.0
## [37] bitops_1.0-6 assertthat_0.2.1 promises_1.1.0
## [40] scales_1.2.1 nnet_7.3-12 gtable_0.3.0
## [43] phylobase_0.8.6 rlang_1.0.1 genefilter_1.68.0
## [46] splines_3.6.2 lazyeval_0.2.2 acepack_1.4.1
## [49] checkmate_2.0.0 yaml_2.2.1 reshape2_1.4.3
## [52] crosstalk_1.0.0 backports_1.1.5 httpuv_1.5.2
## [55] Hmisc_4.3-1 tools_3.6.2 gridBase_0.4-7
## [58] ellipsis_0.3.2 jquerylib_0.1.3 Rcpp_1.0.6
## [61] plyr_1.8.5 base64enc_0.1-3 progress_1.2.2
## [64] zlibbioc_1.32.0 purrr_0.3.3 RCurl_1.98-1.1
## [67] densityClust_0.3 prettyunits_1.1.1 rpart_4.1-15
## [70] pbapply_1.4-2 viridis_0.5.1 ggrepel_0.8.1
## [73] cluster_2.1.0 here_0.1 magrittr_1.5
## [76] data.table_1.12.8 RSpectra_0.16-0 RANN_2.6.1
## [79] hms_0.5.3 mime_0.9 evaluate_0.14
## [82] xtable_1.8-4 XML_3.99-0.3 jpeg_0.1-8.1
## [85] sparsesvd_0.2 HSMMSingleCell_1.6.0 compiler_3.6.2
## [88] tibble_3.1.6 crayon_1.3.4 htmltools_0.5.1.1
## [91] mgcv_1.8-31 later_1.0.0 Formula_1.2-3
## [94] tidyr_1.0.2 howmany_0.3-1 DBI_1.1.0
## [97] MASS_7.3-51.4 Matrix_1.3-2 ade4_1.7-13
## [100] cli_3.1.1 igraph_1.2.4.2 pkgconfig_2.0.3
## [103] bigmemory.sri_0.1.3 rncl_0.8.3 registry_0.5-1
## [106] foreign_0.8-72 locfdr_1.1-8 xml2_1.3.2
## [109] foreach_1.4.7 annotate_1.64.0 bslib_0.2.4
## [112] rngtools_1.5 pkgmaker_0.31 webshot_0.5.2
## [115] XVector_0.26.0 bibtex_0.4.2.2 stringr_1.4.0
## [118] digest_0.6.27 DDRTree_0.1.5 softImpute_1.4
## [121] rmarkdown_2.11 htmlTable_1.13.3 edgeR_3.28.0
## [124] kernlab_0.9-29 shiny_1.4.0 lifecycle_1.0.1
## [127] monocle_2.14.0 nlme_3.1-142 jsonlite_1.6.1
## [130] Rhdf5lib_1.8.0 viridisLite_0.3.0 limma_3.42.1
## [133] fansi_0.4.1 pillar_1.6.4 lattice_0.20-38
## [136] fastmap_1.0.1 httr_1.4.1 survival_3.1-8
## [139] glue_1.4.2 qlcMatrix_0.9.7 FNN_1.1.3
## [142] png_0.1-7 iterators_1.0.12 bit_1.1-15.2
## [145] stringi_1.4.5 sass_0.3.1 HDF5Array_1.14.1
## [148] blob_1.2.1 latticeExtra_0.6-29 memoise_1.1.0
## [151] dplyr_1.0.3 irlba_2.3.3 ape_5.3